%% CTBCS 3D Benettin‑QR
clc; clear; close all

%—— 参数设置 ——————————————————————————————
r1_list   = linspace(-20,20,400);   % 扫描区间
b1        = -3;                     % 移位常数
X0        = [0.3; 0.5; 0.6];        % 初始向量 [x0;y0;z0]
N_trans   = 300;                    % 暂态迭代步数
N_iter    = 500;                    % 后续迭代步数（用于 LE 计算）
deltaDiff = 1e-7;                   % 数值差分扰动步长
%———————————————————————————————————————————————

LEs_all = zeros(numel(r1_list),3);

parfor k = 1:numel(r1_list)
    LEs_all(k,:) = benettinQR_CTBCS3D(r1_list(k), b1, X0, N_trans, N_iter, deltaDiff);
end

%—— 绘图 —————————————————————————————————————
figure
plot(r1_list,LEs_all(:,1),'r','LineWidth',1.2); hold on
plot(r1_list,LEs_all(:,2),'g','LineWidth',1.2);
plot(r1_list,LEs_all(:,3),'b','LineWidth',1.2);
xlabel('r_1'); ylabel('LEs');
xlim([-20 20]); ylim([-10 15]);
set(gca,'FontSize',13)
title('(b)CTBCS LEs (复数)','FontWeight','normal')
legend('LE_1','LE_2','LE_3','Location','Best')
grid on
%———————————————————————————————————————————————

%% ====== 本文件所有子函数，下移到脚本末尾 ======

function LEv = benettinQR_CTBCS3D(r1, b1, X0, N_trans, N_iter, delta)
% 3D CTBCS + Benettin‑QR 最大Lyapunov指数计算
    % 1) 丢弃暂态
    X = X0;
    for i = 1:N_trans
        X = CTBCS3D_step(X, r1, b1);
    end
    % 2) 初始化 QR 基和累计变量
    Q     = eye(3);
    sumLE = zeros(3,1);
    % 3) Benettin‑QR 核心循环
    for i = 1:N_iter
        J = CTBCS3D_jacobian(X, r1, b1, delta);
        Z = J*Q;
        [Q,R] = qr(Z,0);
        sumLE = sumLE + log(abs(diag(R)));
        X = CTBCS3D_step(X, r1, b1);
    end
    % 4) 平均并从大到小排序返回
    LEv = sort(real(sumLE / N_iter), 'descend')';
end

function X2 = CTBCS3D_step(X1, r1, b1)
% CTBCS 三维一次迭代：F₃ 改量子Logistic，G₃为3D-SIMM
    % （1）种子系统 F3：量子 Logistic 3D
    F3 = quantum3D(X1, 3.99, 30);
    % （2）种子系统 G3：3D-SIMM
    G3 = SIMM3D_vec(X1, 1, 2*pi, 11.5);
    % （3）CTBCS 复合映射
    X2 = cos(pi*(r1*F3 + (1-r1)*G3) - b1);
end

function J = CTBCS3D_jacobian(X, r1, b1, delta)
% CTBCS 三维雅可比（链式法则 + 数值差分）
    % 基态与映射值
    F3 = quantum3D(X,3.99,30);
    G3 = SIMM3D_vec(X,1,2*pi,11.5);
    Arg = pi*(r1*F3 + (1-r1)*G3) - b1;     % 3×1
    C   = -sin(Arg);                       % d(cos)/d(arg)
    % 数值差分求两套种子系统雅可比
    JF = jacobian_numeric(@(Y) quantum3D(Y,3.99,30), X, delta);
    JG = jacobian_numeric(@(Y) SIMM3D_vec(Y,1,2*pi,11.5), X, delta);
    % 组合： J = diag(C) * π*(r1*JF + (1-r1)*JG)
    J = diag(C) * (pi*(r1*JF + (1-r1)*JG));
end

function J = jacobian_numeric(f, X, delta)
% 通用数值差分雅可比
    n  = numel(X);
    f0 = f(X);
    J  = zeros(n,n);
    for j = 1:n
        Xp = X; Xp(j)=Xp(j)+delta;
        f1 = f(Xp);
        J(:,j) = (f1 - f0) / delta;
    end
end

function Y = quantum3D(X, r, b)
% 量子 Logistic 三维系统（论文 Eq.(1) 原式实现）
% X = [x; y; z]，r, b 为系统参数
    x = X(1); y = X(2); z = X(3);
    % 复数共轭 / 点上标表示
    xd = conj(x);    % 即 xi̇
    zd = conj(z);    % 即 zi̇

    Y = zeros(3,1);
    % xi+1 = r*(xi - |xi|^2) - r*yi
    Y(1) = r*( x - abs(x)^2 ) - r*y;

    % yi+1 = -yi e^{-2b} + e^{-b} r [ (2 - xi - xi̇) yi  -  xi zi̇  -  xi̇ zi ]
    Y(2) = -y * exp(-2*b) + exp(-b)*r*( (2 - x - xd)*y  -  x*zd  -  xd*z );

    % zi+1 = -zi e^{-2b} + e^{-b} r [ 2(1 - xi̇ zi) zi  -  2 xi yi  -  xi ]
    Y(3) = -z * exp(-2*b) + exp(-b)*r*( 2*(1 - xd*z)*z  -  2*x*y  -  x );
end


function Y = SIMM3D_vec(X, a, b2, c)
% 3D‑SIMM 系统（论文 Eq.(2)）
    x = X(1); y = X(2); z = X(3);
    x1 = sin(b2*z)   * sin(c/x);
    y1 = sin(b2*x1) * sin(c/y);
    z1 = sin(b2*y1) * sin(c/z);
    Y  = a * [x1; y1; z1];
end

% 
% %%  时间差分版本
% clc; clear; close all
% 
% %—— 参数扫描 ——————————————————————————————————
% r1_list   = linspace(-20, 20, 400);
% b1        = -3;
% X0        = [0.3; 0.5; 0.6];
% N_trans   = 300;
% N_iter    = 500;
% deltaDiff = 1e-7;
% %———————————————————————————————————————————————————
% 
% LEs_all = zeros(numel(r1_list), 3);
% 
% parfor k = 1:numel(r1_list)
%     LEs_all(k,:) = benettinQR_CTBCS3D_diff(r1_list(k), b1, X0, N_trans, N_iter, deltaDiff);
% end
% 
% %—— 绘图 —————————————————————————————————————
% figure
% plot(r1_list,LEs_all(:,1),'r','LineWidth',1.2); hold on
% plot(r1_list,LEs_all(:,2),'g','LineWidth',1.2);
% plot(r1_list,LEs_all(:,3),'b','LineWidth',1.2);
% xlabel('r_1'); ylabel('LEs');
% xlim([-20 20]); ylim([-10 15]);
% set(gca,'FontSize',13)
% title('(b) CTBCS LEs (时间差分)','FontWeight','normal')
% legend('LE_1','LE_2','LE_3','Location','Best')
% grid on
% %———————————————————————————————————————————————————
% 
% %% ———— 子函数 ————
% 
% function LEv = benettinQR_CTBCS3D_diff(r1, b1, X0, N_trans, N_iter, delta)
% % Benettin‑QR with time‑diff F3 and standard G3
%     % 丢弃暂态
%     X = X0;
%     for i = 1:N_trans
%         X = CTBCS3D_step_diff(X, X, r1, b1);  % 初始两步都用 X_prev = X0
%     end
% 
%     Q     = eye(3);
%     sumLE = zeros(3,1);
%     for i = 1:N_iter
%         % 解析数值 Jacobian
%         JF = jacobian_numeric(@(Y) quantum3D_diff(Y,Y,3.99,30), X, delta);
%         JG = jacobian_numeric(@(Y) SIMM3D_vec(Y,1,2*pi,11.5),  X, delta);
%         % 复合 CTBCS Jacobian
%         F3  = quantum3D_diff(X, X, 3.99, 30);
%         G3  = SIMM3D_vec(X, 1, 2*pi, 11.5);
%         Arg = pi*(r1*F3 + (1-r1)*G3) - b1;
%         C   = -sin(Arg);
%         J   = diag(C) * (pi*(r1*JF + (1-r1)*JG));
% 
%         [Q,R] = qr(J*Q,0);
%         sumLE = sumLE + log(abs(diag(R)));
% 
%         % 迭代
%         X = CTBCS3D_step_diff(X, X, r1, b1);
%     end
% 
%     LEv = sort(real(sumLE / N_iter), 'descend')';
% end
% 
% function X2 = CTBCS3D_step_diff(X_cur, X_prev, r1, b1)
% % CTBCS 三维迭代，F3 用时间差分，G3 用标准 3D‑SIMM
%     F3 = quantum3D_diff(X_cur, X_prev, 3.99, 30);
%     G3 = SIMM3D_vec(X_cur,      1, 2*pi, 11.5);
%     X2 = cos(pi*(r1*F3 + (1-r1)*G3) - b1);
% end
% 
% function J = jacobian_numeric(f, X, delta)
%     n  = numel(X);
%     f0 = f(X);
%     J  = zeros(n,n);
%     for j = 1:n
%         Xp = X; Xp(j)=Xp(j)+delta;
%         f1 = f(Xp);
%         J(:,j) = (f1 - f0)/delta;
%     end
% end
% 
% function Y = quantum3D_diff(X, X_prev, r, b)
% % 三维量子Logistic，˙x_i 用时间差分
%     x  = X(1); y  = X(2); z  = X(3);
%     x0 = X_prev(1); z0 = X_prev(3);
%     dx = x - x0; dz = z - z0;
%     Y = zeros(3,1);
%     Y(1) = r*(x - x^2) - r*y;
%     Y(2) = -y*exp(-2*b) + exp(-b)*r*((2 - x - dx)*y - x*dz - dx*z);
%     Y(3) = -z*exp(-2*b) + exp(-b)*r*(2*(1 - dx*z)*z - 2*x*y - x);
% end
% 
% function Y = SIMM3D_vec(X, a, b2, c)
% % 3D‑SIMM 原式
%     x = X(1); y = X(2); z = X(3);
%     x1 = sin(b2*z)   * sin(c/x);
%     y1 = sin(b2*x1) * sin(c/y);
%     z1 = sin(b2*y1) * sin(c/z);
%     Y  = a*[x1; y1; z1];
% end
